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QQ ' Abstract 

^^ ' We present an algorithm to solve — Au — f{x, u) = g with Dirichlet bound- 

ary conditions in a bounded domain Q.. The nonlinearities are non-resonant 

Q, ■ and have finite spectral interaction: no eigenvalue of —Ad is an endpoint of 

,^f^ ' 92/(f2,R), which in turn only contains a finite number of eigenvalues. The 

. , algorithm is based in ideas used by Berger and Podolak to provide a geomet- 

,S^ • ric proof of the Ambrosetti- Prodi theorem and advances work by Smiley and 

jrt \ Chun for the same problem. 
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00 ■ 1 Introduction 

; . In this paper, we describe a numerical algorithm to solve 

O ; - ^u{x) - /(.T, u{x)) = 5(x), u|ao = 0. (1) 

Here, the nonlinearity / is an appropriate function, to be defined later, and the 
domain i7 G K" is open, bounded, connected and has Lipschitz boundary dVl. 
\^u ' Theoretical tools go hand in hand with numerical methods. Local behavior at 

'V^ . regular points concerns both the inverse function theorem and Newton's inversion 

C^ ' algorithm. Homotopy arguments such as degree theory go along well with continu- 

ation methods. The celebrated mountain pass lemma f|15|) is the starting point of 
an algorithm presented in [B] . More recently, ideas used in computer assisted proofs 
were combined with the topological toolbox with striking effect by Breuer, McKenna 
and Plum ([5]). In this paper, we explore numerically a global Lyapunov-Schmidt 
decomposition. 

More precisely, we consider a class of C^ maps F : X ^ Y between Banach 
spaces. Split X = Wx ® Vx and Y = Wy © Vy into closed horizontal and vertical 
subspaces. Define complementary projections Py,Qy : Y ^ Y so that Ran Py = 
Wy and Ran Qy = Vy . A map F is flat if, for each x e X, Py o F : x + Wx -^ Wy 
is a diffeomorphism. Here x + Wx is the affine horizontal subspace obtained by 
translating Wx by x. 

This stringent hypothesis gives rise to substantial geometric structure. Images 
under F of affine horizontal subspaces, sheets, intercept transversally affine vertical 
subspaces y + VY,y € Y, ed, a unique point. Preimages under F of affine vertical 
subspaces, fibers, are submanifolds of X diffeomorphic to Vx- Indeed, X is foliated 
by fibers, and each fiber intersects transversally each affine horizontal subspace at 
a unique point. 



The bifurcation equations related to the decomposition for F{x) = y are 

P^F{w + v) = P,.y, Q,.F{w + v) = Q^y, weW^, v e Vj, 

and the first equation, by flatness, admits a unique solution w{v) for each fixed v. 
Said differently, given y G Y, w{v) is the unique point of the fiber F^^{y + Vy) 
in the affine horizontal space v + Wx- Clearly, the fiber through w(w) contains all 
solutions of F{x) = y. 

In a nutshell, the algorithm first computes w{v) based on the finite element 
method. The search for solutions of F{x) = y then reduces to inverting a (com- 
putable) map between isomorphic finite dimensional subspaces Vx and Vy 

We use piecewise linear finite elements. For the sake of sparsity, we exploit 
the decompositions of spaces X and Y. Indeed, from flatness, a large part of the 
derivative DF{u), taking Wx to WV, is invcrtible. We extend this isomorphism to 
one from X to Y which is especially simple to code. The search for w{v) becomes 
then a standard continuation method between given affine horizontal subspaces, 
with the advantage that computations are performed in the full spaces X and Y . 

The restriction of -F to a fixed fiber a now can be computed by a predictor- 
corrector algorithm. In more detail, Vx parametrizes a and, given x S a, a point 
X + v,v G Vx corresponds to a unique point x € a with the same height (i.e., 
Qxix + v) — Qxx), obtained from x -f i; by the same continuation method used to 
compute w{v). We are then left with inverting a (constructible) function between 
finite-dimensional spaces. 

The constructions rely on the assumption that the projections Qx and Qy are 
computable. In order to avoid unnecessary abstraction, we present the algorithm 
for the special case related to equation ([!]). This allows us to discuss some imple- 
mentation issues. We consider the map F{u){x) = —Au{x) — f{x,u{x)) between 
Sobolev spaces X = H^{^) and Y = H-^{n) ~ -ffd(^). The nonhnearity / is 
assumed to be appropriate: 

feC^nxR), ||/(-, 0)11^2 <oo and ||a2/|U < oo. 

Here 82 is the partial derivative with respect to the second variable. Let / — [a, b] be 
an interval containing 82/ {^,M.). Set Vx — Vy to be the direct sum of the (maximal) 
invariant subspaces associated to the eigenvalues in [a, 6] of — A^i. Finally, let 
Wx = V^ and Wy = Vj- in X and Y respectively. As we shall see in Section [3l 
F : X — > F is flat with respect to this decomposition. 

The search for w{v) is robust and globally stable: errors self-correct in the spirit 
of Newton-type iterations, and the linear operators which require inversion are both 
uniformly bounded and uniformly coercive. Searching for solutions in the fiber is 
not necessarily an easy task. When dimV^ = dimVV = I7 one needs to invert a 
function from R to R, and root solvers abound. For higher dimensions, matters are 
harder. For the two-dimensional case, we present an example in Section [6.31 

The history of semilinear elliptic theory, together with some computational as- 
pects, is very well described in [5]. Here we emphasize some techniques which are 
relevant to our text. A good introduction to computer assisted proofs in a related 
context is [13] ■ 

Hammerstein ([5]) and Dolph CS]) showed that, if 92/(il,R) does not contain any 
eigenvalue of —Ajj, the map F : X —i' Y is a. (global) diffeomorphism. For the choice 
Vx = Vy = {0}, F is trivially flat. Numerical inversion might proceed by standard 
continuation methods, requiring inversion of —Av{x) — d2f{x,u{x))v{x) = h{x)^ 
with Dirichlet boundary conditions. 

In the autonomous case f(x,u) = f{u), Ambrosetti and Prodi ([I,) presented a 
thorough analysis of a situation in which F{u) = g admits multiple solutions. Their 
result, immediately amplified by Manes and Micheletti ([12]), essentially states that 



if / is convex and /'(K) only contains the smallest eigenvalue Ai of —Ad, then 
g can only have 0, 1 or 2 preimages. Later, the Ambrosetti-Prodi map F was 
given a novel geometric description by Berger and Podolak ([4 ): for the Lyapunov- 
Schmidt decomposition Vx — Vy — Span{(y9i}, they showed that F is flat. Here 
Lpi is a positive eigenfunction associated to Ai. Their proof uses a coercive bound 
\\DF{w + w)|| > C||it;||, uniform in v. They then showed that each fiber a, the 
preimage under F of a vertical affine subspace, is a differentiable curve. Moreover, 
the restriction of F to each a becomes s i— )► — s^, after a change of variable. Since 
fibers foliate X, F is a global fold: (global) changes of variables from X and Y 
to a common space Z convert F into F : Z ^ Z where Z = V ®W given by 
F{s,r) i-> {-s'^,r). 

Hess ([10]) extended the result by Ambrosetti and Prodi to the nonautonomous 
case. The Lyapunov-Schmidt decomposition in this context seems to have been es- 
tablished by Smiley and Chun, who also realized its potential for numerics: solving 
F{u) = g boils down to solving the equation restricted to each fiber ([H], [T7].|18|). 
In these papers, the authors are concerned with approximating the bifurcation equa- 
tions, i.e., the restriction of F{u) = g to the fiber ag whose image contains g, using 
finite element methods. To solve the inversion problem of F restricted to a given 
fiber. Smiley and Chun ([12]) developed a general solver for locally Lipschitz maps 
from K" to R", and provided examples ([20])- Our algorithm, on the other hand, 
computes a point in ag, for arbitrary g and fixed affine horizontal subspace. As a 
byproduct, it yields a (stable) procedure to move along ag. 

In [13], Podolak considered fibers for different nonlinearities. In jll], fibers were 
used to show that the map G{u{t)) = u'{t) + v?{t) — u{t) is a global cusp from 
the space of periodic functions in C^([0, 1]) to C([0, 1]): after global changes of 
variables, G becomes (x, y) i— ;> (x^ — xy, y). 

In Section [21 the consequences of flatness are presented in the general setting 
of a map between Banach spaces. The techniques are standard and the reader is 
invited to skip the section if he feels comfortable with the implications of flatness 
stated above. For equation ([Ij, flatness follows from the coercive bound proved in 
Section [3] The algorithm is described, first theoretically and then in more concrete 
terms, in Sections |4] and [5] We finish with some examples in Section (6] 

The authors gratefully acknowledge support from CAPES, CNPq and FAPERJ. 

2 Geometry of flat maps 

Let X and Y be Banach spaces which split as direct sums of horizontal and vertical 
subspaces, X — Wx © Vx and Y = Wy ® VV- We assume all four subspaces to be 
closed and define the pairs of (bounded) complementary projections Px + Qx = Ix 
and Py + Qy = Ix, where Px and Py (resp. Qx and Qy) project on horizontal 
(resp. vertical) subspaces. Sets of the form x + Wx (resp. y + VFy) or x + Vx (resp. 
y + Vy) will be denoted by horizontal and vertical affine subspaces. 
For V GVx, the projected restriction 



F,:Wx^Wy, F, (w) = PyF{w + v) 



acts between horizontal subspaces. A C-^ map F : X —^ Y is flat with respect to 
a decomposition of X and Y as above if, for any v & Vx, the associated projected 
restriction F^ is a diffeomorphism. Thus, F takes horizontal affine subspaces x + Wx 
injectively to their images, which are graphs of functions from Wy to Vy'. the 
surfaces F{x + Wx) are called sheets. 

The situation is familiar: horizontal variables are trivialized by a change of 
variables and vertical variables are the unknowns of the bifurcation equations. This 
is the content of the next proposition. 



Proposition 1 Let F : X ^ Y be flat for decompositions of X and Y as above. 
Then the function 

<^:X^W^®V^^W^(BV,, $(z,«) = ((F,)-i(z),«) 

is a C^ difjeomorphism such that F = Fo^: X^s-Y is F(z, v) = (z, 0(z, v)) for a 
C^ function (j) : X ^f Vy ■ 

Proof: As usual in Lyapunov-Schmidt arguments, consider ttie tentative inversion 
of an arbitrary point {wy, Vy) ^ Wy © Vy, 

F{w + v) = {Fy{w), QyF(w,v)) ^{wy,Vy). 

By hypothesis, w = (F„)-i(wy). Clearly, $ = {{Fy)-\id) : Wy ®V^ ^ W^ ®Vx 
and ^ = ^~^ are C^ diffeomorphisms. The rest follows from the diagram below. 

F 

(w. v) I y (wy : Vy) 



{wy : V) 

u 

A fiber a is the preimage of a vertical afSne subspace. We denote by ag the 
fiber which is the preimage of the affine subspace Vy + g. The height of a point 
X E X (rcsp. y E Y) is the vector QxX (resp. QyTj)- 

Proposition 2 Let F : X ^ Y be flat. Then each fiber ag is a C^ surface of 
dimension dim Vx, which intersects each horizontal affine subspace at a unique point 
X transversally, i.e., X = Txag(BWx ■ The height map x i-> QxX is a diffeomorphism 
between the fiber ag and the vertical subspace Vx, with inverse T-Lg :Vx ^ Q.g given 
byngiv)=v + F-^PYg. 

According to the proposition, Wx parametrizes (bijectively) the set of fibers, 
and Vx each fiber. Horizontal affine subspaces are sent injectively by F to their 
images but fibers are not necessarily taken injectively (nor surjectively!) to ver- 
tical subspaces. In particular, the given hypotheses are not enough to imply the 
properness of the map F : X ^>^ Y. 

Proof: We use the notation from the previous proposition. The change of variables 
<i>(z, v) = ((_F'u)~^(z), v) is a diffeomorphism from each vertical affine subspace in X 
to a fiber of F with the property that heights are preserved. Each statement about 
fibers follows easily from its counterpart for vertical affine subspaces in X. ■ 

Proposition 3 Let F : X ^)- Y be flat. Then sheets are manifolds which intersect 
vertical affine subspaces transversally. If Xc is a critical point of F contained in the 
fiber a, then KeT{DF{xc)) C Tx^a. 

Transversal intersection at F{x) e Y means Y — DF{x)Wx ® Vy- 
Proof: The change of variables <I> is a diffeomorphism between horizontal affine 
subspaces, so sheets of F are also the images of horizontal affine subspaces Wy 
under F, and hence are manifolds of codimension dimW. Clearly, the tangent 
space of a sheet at a point F{x) consists of the closed vector space DF{x)Wx. To 
see that DF{x)Wx CiVy — {0}, suppose DF[w + v)w = v,ioT x = w + v and w,w E 
Wx,v eVx,v e Vy. Then PyDF{w + v)w = and since PyDF{w + v) = DFy{w), 
we have that DFy(w)'w = 0. Now, F^ is a diffeomorphism between v + Wx and WV 
and thus w = 0. Counting dimensions, we conclude that Y = DF(x)Wx © Vy. 

At a critical point Xc E a, Py o DF{xc) '■ Wx -^ Wy is an isomorphism by 
ffatness, and thus DF{xc) ■ Wx -^ DF{xc)Wx also is. Split X = Wx ® T^^a as in 
Proposition^ An element of the kernel of DF{xc) : Wx ^T^^a ^ DF{xc)Wx ffl Vy 
must have null cooordinate in Wx, so that KeT{DF{xc)) C Tx^a. ■ 



Thus, the projection Py is a diffeomorphism between each sheet and Wy- Fibers 
are disjoint, but sheets are not — this is why some points have more than one 
preimage under F. 

3 Smoothness and Flatness 

Let Q, denote an open, connected, bounded set of M" with Lipschitz boundary 
dVt. We use standard notation for Sobolev spaces W^''p{VI) , Wg'^(ri); i.e., the j- 
seminorm of a function u G W'''P{fl) is given by |m|^.p = J2\a\=i i-^""llp' ^^"^ i*^ 
j-norm by ||it||^p = X]fc<i I'^lfc.p- Of special interest are the spaces with p = 2 and 
j = 0, 1 or 2, which we~denote by H°{n) = L'^{n) , H^{n) and H^{n). In the case 
of Dirichlet boundary conditions we work mainly with the space HQ^fl) or with 
i/o(^)- Using Poincare's inequality we see that seminorms of the H spaces are 
indeed norms, equivalent to the full norms, and that they are Hilbert spaces, with 
inner products given by 



{u,v)o— / uv, {u,v)i ^ / Vu ■ Vw, {u,v)2— / Au-Av. 
Jn Jn Jn 

We identify i?o(ri) ~ H~^{n) via {u, •) = {u, •)i, where the tilde denotes the func- 
tional induced by an element oi Hg(fl) and the brackets (, ) with no subscript denote 
the coupling between a space and its dual. Thus 

(m, «)_! = (u, f )i, ||u||_i = |w|i, u„ — > u -i^ Uji —^ u and u i— > (m,-)i. 

Recall that a function /ifixK— >]Risa Caratheodory function if s h-> /(a;, s) is 
continuous for almost every x and x ^-^ f{x, s) is measurable for every s. A map of 
the form u i-t- f{-,u) is a Nemytskii operator. We work with appropriate functions 
/, which satisfy 

feC\nxR), ||/(.,0)||o<(X) and f^a/IU < oo. 

Set X = iJQ(ri), Y = H^^{fl) and, for an appropriate /, the nonlinear map 

F:X~^Y, F{u){x) ^ -Au{x)- f{x,u{x)). 

The (Dirichlet) Laplacian acts weakly and /(-,«(•)) is the functional given by 

z^ {fi-,u),z)o^ f{x,u{x))z{x)dx. 



Write F{u) = —Au — Nf{u), where Nf{u){x) = f{x,u{x)) is the Nemytskii 
operator associated with /. The estimates below are a minor extension of the 
properties enumerated in [5]. Throughout the text, C denotes a positive constant, 
which may change along the argument. 

Proposition 4 Let f be an appropriate function. Then F : X -^ Y is a C^ map. 

Proof: It suffices to show that Nf : H^{il) — s- L'^(il) is C^. Indeed, this implies 
that Nf : X -^ Y is also C\ since X C H^{n) and L'^{fl) C Y. Notice first that 
Nf : H^{il) -^ L^(ri) is well defined. Indeed, by the Taylor formula with integral 
remainder in x there exist a, 6 > with 

\\fi;u)l<\\al + b\\ul<Cil + \\ul). 



Write the superlinear remainder 
e(x, h{x)) — f{x, u{x) + h{x)) — f{x, u{x)) — d2f{x, u{x))h{x) ~ 5(x, h(x)) h(x) 

where 

5{x,h{x)) :— i d2f{x,u{x)+Th{x))~d2f{x,u{x)) dr. 
Jo 

To ensure differentiabihty, we need to show that ||e||o ^- as \\h\\i -^ 0. By the 
Sobolev imbedding theorems, h G -ff^(il) is also in L*(ri) for some s > 2 (if n > 2 
we can take any 2 < s < 2* = -^^: whereas for n < 2 any s > 2 works). By Holder's 
inequality, ||e||o < ||/i||o,3||<5||o,r, where r > 2 is given by ^ + ^ = i. Again from the 
imbedding theorems, ||/i||o,3 < C'l^lu and we are left with showing that ||5||o.r -^ 
as ||ft-||i -^ 0. Switching to a subsequence if necessary, h —> pointwise a.e. so 
that the integrand |5|^'' also converges to zero pointwise a.e., by the continuity of /. 
Hence, by the bounded convergence theorem, ||(5||o,r — >■ 0. This holds for any sub- 
subsequence and thus for \\h\\i -^ in general. That z i-t- d2f{-,u)z is a bounded 
map follows from \\d2f{-,u)z\\n < ||92/|U||2;||o < C||(92/|U||2;||i, completing the proof 
of Frechet-differentiability. 

We now show continuity of the derivative: 

For any u e iJi(f^), \\DNf{u + h)-DNf{u)\\^0 whenever \\h\\, ^ 0. 

Suppose V e H^{Q), set g{x,h{x)) — d2f{x,u{x) + h{x)) — d2f{x,u{x)) and notice 
that, for r, s defined above, 

\\9{;h)vl<\\gi-.h)\U\vl,.<C\\g{;h)lJvl. 

The constant C does not depend on the function v and, since / is C^, the previous 
argument with the bounded convergence theorem holds, yielding 

\\Nf{u + h)-Nf{u)\\^0, as||/i||i^O. 



In Theorem [I] we exhibit direct sums X = Wx ® Vx and Y = WV © V^ for 
which F : X — > y is flat. We follow closely some arguments in \T.. 

Denote the (possibly repeated) eigenvalues of — A/j : -ffo (i^) C H'^{fl) — > iy°(0) 
by < Ai < A2 < . . . and choose corresponding orthogonal cigenfunctions ipk,k = 
1,2,.. .. Orthogonality holds for the four spaces iJg (17), H^in), H°{n) and H-^{n). 
Now let / = [a, b] an interval containing 82/ {^, M). The I -index set is I = {i \ \i E 
I}. The I -decompositions of X and Y are defined as follows. Set Vx — Vy be the 
subspace spanned by {(p^, i G I}. Also, set Wx — Vx and Wy — Vy ■ The four 
projections Px,Qx,Py and Qy are now orthogonal. As in the previous section, for 
V G Vx, the projected restriction F-„ : Wx -^ Wy acts between horizontal subspaces. 



Proposition 5 Let f be an appropriate function, I = [a, 6] D 92/(r2,R) and I- 
decompositions X = Wx © Vx and Y — Wy © Vy ■ Then the derivatives DFy of the 
associated restricted projections of F : Wx © Vx — >■ Wy © Vy are uniformly bounded 
from below. More precisely, there exists C > such that 

Vw e Vx Vw e Wx V/i G Wx, \\DF^{w)h\\_, > C\\hl. 

Also, all derivatives DFy are invertible. 

When / — [a, b] contains only the first eigenvalue Ai, i.e., 2 — {1}, this estimate 
has been extensively used ([1], [4]). It is also used in [7] in the case when / contains 



the first k eigenvalues of — A/^. Tlie nonautononious case has been considered in 
[TU] and [12]. The result below is slightly more general. 

Proof: From Proposition 01 each restricted projection Fy : Wx -^ WV is C^ with 
derivative Z? Ft, (w) : Wx — > Wy givenhy DFy{w)h{x) — -Ah{x)—PYd2f{x,u{x))h{x), 
where u = w + v. Take h e Wx of unit norm and let "f — (a + b)/2. Adding and 
subtracting 7/1, 

\\DFyiw)h\\_, = \\PYi-Ah - -fh) - P^(92/(-, u)h - 7/i)||_i 

> iP,(-A/i - 7/i)||_i - ||P^(a2/(-, n)h - ^h)\_, 
>\Ah\_,-\Bh\_,. (2) 

We bound ||_B/i||_i from above. For z <^ X and w G Wx-, 

\Bh\_^ =sup {PY{d2f{-,u)h-'-fh),z) = sup {d2f{-,u)h-jh,w) 

\\z\\ = l ||»|| = i 

= sup {{d2f{-,u) ~j)h,w)„ < \\d2f{-,u) -7||<^ sup {\h\,\w\)o. 

\\w\\ = l \\w\\ = l 

By Cauchy-Schwartz, the supremum is realized when \w\ is a scalar multiple of \h\, 
which is the case when w = p/i, p 6 K. Since w and h are unit vectors in X, we 
may take p = 1 and, defining c = ||52/(-, u) — j\\^, 

||i?/i|U. <c(|/l|,|/.|)„ = c||/.||2 = ^c/^^||^,L^ = ^(c/A,)/i^||^,||?, (3) 

where hk G K are the coefficients of the expansion of h in eigenfunctions and I — 
{£<...< r} is the /-index set. To estimate ||A/i||_i from below, start with 

\\Ah\\_i = sup (Py-(-A/i-7/i),z) = sup (-A/i -7/1,10) 
lkll=i ll^lhi 

= sup {{h,w)i -^{h,w)o). 

Split VFx = W^_ © W+, where 

W- = {u : u = '^ Ukifik}, W+ ^ {u : u = '^ Ukifik} 

k<l k>r 

are orthogonal subspaces in the H^ and i/° norms. Split h = h^ + h+ and set 
w = h^ — h- G X, clearly a unit vector: 

IIA/^IU, > {h,h+ - h_), - ^{h,h+ - h_l = (|/^+|2 - 7||;i+||2) + (^||/j_||2 _ |/,_|2) 

fc>r fc<i 

fc>r fc<i 

Notice that (1 — •y/Xk) is positive (resp. negative) for /c > r (resp. k < i). Then 



E(i - i/^k)hi\^k\', + Y^h/Xk i)hi\^k\i 



\\AhU > E 11 - ^/^>'\ '^'l^fcl? = T.(^k/X,)hl\^k\l (4) 

where Ck = \Xk — 7|. Combining equations ([2|), ([3]) and (|4|), 

\\DF^{w)h\\_, > E(Cfe - c)/A, /i2|^,|2 > [myCfe - c)/Afc] E /^fcl^fel? 

inf^(Cfe - c)/Afe') |/j|? = C|/i|? = C. 



The infimum above is achieved at one of the outer eigenvalues closest to [a, b], prov- 
ing the injectivity of DFy{w). We now show that DFy{w) is a Fredholm operator 
of index zero, and hence surjective. 

Indeed, —A : X -^ Y is an isomorphism and DF{u) : X ^ Y given by 
DF{u)z = —Az — d2f{-, u)z is obtained by adding a compact operator, from Propo- 
sition HI Now, DFy{w) — Py o DF{v + w) o l, where the projection P^ : F — ;■ W^ 
and the inclusion t : Wx — >■ X are Fredholm operators, whose indices add to zero. 
Thus DFy(w) is also Fredholm of index zero. ■ 

Recall Hadamard's global inversion theorem ([3]). 

Lemma 1 Let $ : X — > y be a C^ map between Banach spaces X and Y such that 
D^{u) is invertible for each u Cz X . Suppose there exists C > such that 

yu,hex ||i:i$(u)/i|| >c||/i||. 

Then ^ is a global C^ -diffeomorphism. 

Theorem 1 Let f be an appropriate function, I — [a,b] D 92/(^,R)- Then the 
map F : X ^f Y is flat with respect to the I -decompositions X — Wx ® Vx and 

Y ^Wy®Vy. 

Proof: Simply combine the proposition and lemma above. B 

There is an analogous statement for F : Hq — ;> H^. 

From the previous section, since F is flat, its domain is foliated by C^ fibers of 
dimension dim Vx = dim Vy , which are transversal to the horizontal affine subspaces 
X + Wx and are parameterized diffeomorphically by the height function. The bound 
in Proposition [S] allows to make precise the idea that fibers are uniformly steep and 
sheets are uniformly flat. 



Proposition 6 Let f be appropriate, I = [a, 6] D d2f{i^,R) and Wx ® Vx and 
Wy ffi Vy be the corresponding L -decompositions of X and Y . Denote by X an index 
set for the basis of eigenf unctions spanning Vx ~ Vy and let 

u{t)=w(t)+v{t), forw(t)eWx, v(t)^Y.^^^^e^^^ i = (ii,...,t|i|) 

be a parametrization of a fiber a of the flat map F : X ^^ Y . Then there exists a 
constant C , independent of t, such that 

\Vtw{t)l<cY,\\^.,l. 

In particular, there exist constants A, B, independent oft, such that 

\\w{t)\l < A + B||t||. 

Let u G X and consider the sheet VV« — F{u + Wx) with tangent space Tp(^^-)Wu cit 
F(u). Then the angle between a vector in T^(ti)W'u and its orthogonal projection in 
Wy is (uniformly) bounded above by a constant less then 7r/2. 

This result is a source of robustness for the numerics in the next sections. 

Proof: Fibers are inverses under F of vertical afhne subspaces in Y. Taking 
derivatives of Py,F(u(t)) = const, with respect to ti, 

DPYF{u{t)) dtMt) = PYDF{u{t)) dtMt) = PYDF{u{t)) {dtMt) + f^) = 0. 



Since for any h e Wx we have PYDF{u{t))h — DFy(i-f{w{t))h, for h = dt^w{t), 

DF,^t){w{t))dtMt) = PYDF{u{t))dtMt) = -PvDF{u{t))^^. 
Using first the lower bound in Proposition [5] and then the boundedness of DF, 

C\\\dtMt)\\i < \\DF,^t)iw{t))dtMt)\\-i = \\PYDF{u{t))^,\\_, < C2y^i, 

for some constant C2. Thus ||VtK;(i)||i < CJ^tei Wfih^ fo^ some other constant C. 
A bound of the form ||w(i)||i < A + B\\t\\ is now immediate. 

To see that Tpf^jWu is bounded away from the vertical subspace, consider the 
sequence of simple estimates, for h G Wx- 

Ci\\hl < \\PYDF{u)h\\_, < \\DF{u)hU < CsWhl 

The cosine between a vector DF{u)h £ Tp^^-^Wu and the horizontal subspace WV 
is \PYDF{u)h\_^/ \DF{u)h\_i, which is bounded from below by C1/C3. ■ 

A regularity theorem for F, in the sense that g G C°"{^) => u € C°°{il) for 
F{u) = g, would imply that points in the same fiber have the same differentiability. 

4 Finding Preimages under F 

We now describe an algorithm to solve F{u) — — Au — /(•, u) == g, ujgo = 0. The 
details of implementation are handled in Section [S] The equation is interpreted as 
the computation of the preimages of g under F : X = Hq^H.) — )> y = H~^{il.). 

For g G H^, there is an alternative point of view, which we do not treat in 
this paper: one might work instead with F : Hq{Q,) -^ i?"(i7), which shares the 
same geometric properties than F, as commented below Theorem [TJ However, the 
discretizations will be performed by choosing appropriate finite elements, and the 
programming becomes easier for the less restrictive basis used in H^, as opposed to 
the finer elements in H^. Clearly, the preimages of 5 G H^ under F are in H^. 

We assume that the nonlinearity / is an appropriate function and the interval 
/ — [a, b] contains (?2/(i^,R), so that, by Theorem[Tl F : X ^ Y is flat with respect 
to the /-decompositions X = Wx ®Vx and Y ^Wy ®Vy. 

In a nutshell, split g = PYg + QYg — gw + gv- The inversion under F of the 
vertical affine space gw + Vy gives rise to a fiber ag which contains all the solutions 
of the original equation. The algorithm first identifies, for a fixed u G Vx, a point 
Ug € ag n {v + Wx}- this is essentially handling the equation PYF{ug) = gw in 
{w + Wx}. The search for solutions then boils down to a finite dimensional problem 
along ttg, which corresponds to the bifurcation equation QyF{u) — gy- 

4.1 Moving in the space of fibers 

Our first goal is to reach a point Ug in the fiber ag = F^^{g + Vy), or more 
realistically, close to it. For an arbitrary u G Vx, we search the unique point u of 
ag in the horizontal affine space v + Wx given by Proposition [5] This is equivalent 
to solving 

PyF{v + w)= Fy{w) = PYg, weWx- 

Since F is flat, for each v & Vx, F^ : Wx — > WV is a diffeomorphism so that, for any 
w G Wx, DFy(w) = — PyA — P5^92/(-, u) is an isomorphism. Thus, we may consider 
Newton's method on F^ to move horizontally in Wx- However, if we restrict our 
computations to horizontal subspaces, our finite elements discretizations will not 
yield sparse matrices. It is natural, then, to search for an extension to the full space 
of the operator DFy which is invertible and easy to compute. 



For w e X we define the linear operator Lc{u) : X ^ Y hy 
Lc{u)z = -Az - PYd2f{:,u)PxZ - cQyQxZ = -Az - PYd2f{-,u)PxZ - cQxZ, 
since Vx —Vy. 

Proposition 7 Write u = w + v € Wx © Vx- The restrictions of Lc{u) : X ^^ Y 
to Wx and Vx are DF^{w) : Wx ^ Wy and -A - cI : Vx ^ Vy- 

Proof: For z = zw + zy e Wx (SVx, 

Lc{u)z = ~A{zw + zv) - PYd2f{-,u)Px{zw + zy) ~cQx{zw + zv) 
= {~Azw - PYd2f{-,u)zw) + {-Azv -czy) 
= DFy{w)zw + (-A - c/)zy. 



Notice that Lc{u) is an integro-difFerential operator. This is not a problem 
for the finite elements discretization and has an added bonus the preservation of 
sparsity of the relevant matrices. 

To reach Ug G v + Wx in the fiber ag , start with an arbitrary uq £ v + Wx ■ To 
update Un G V + Wx , solve 

Lc{un)h = Pyig - F(un)) and set m„+i = m„ + Px ^• 

The projection in the formula for Un+i is redundant, but it removes possible nu- 
merical errors that might give rise to a nontrivial vertical component when solving 
for h. Actually, from Proposition [71 

Un+i^Un + Pxh, where Lc{un)h ^ g - F{un). (5) 

Numerical errors self-correct, in the spirit of Newton's method: termination 
occurs once the norm of the error e„ = Pyi^g — F(un)) is sufficiently small, yielding 
a point Ug essentially in Ug. 

In principle, convergence is not expected and might require prudence: inversion 
of points along the horizontal segment joining PyF(uq) to Pi-g. This always works 
in exact arithmetic, since F^^ is a C^ diffeomorphism. 

The algorithm above also implements the diffeomorphism "Hg : Vx ^>- otg intro- 
duced in Proposition [2] — it suffices to start from v € Vx and move horizontally 
until Newton's iteration reaches ag. 

4.2 Moving Along a Fiber 

Once Ug € Ug is identified, the original problem reduces to a finite-dimensional 
issue. Said differently, we should invert the restriction of _F to a^, which amounts 
to inverting Jg : Vx — >■ Vy given by J-'g = Qy o F o %g . 

Actually, when a value 'Hg{vo) has been computed, we may compute Hgivo+p), 
for p G Vx, by starting from Hgivo) +p instead oi vq +p and then moving horizon- 
tally with Newton's method until we reach ag. The advantage lies in the fact that, 
for small p, there will be less horizontal displacement with the new initial condi- 
tion. The resulting algorithm is essentially a predictor-corrector scheme. Figure [1] 
illustrates the procedure in the case when dim Vx = |I| = 1, so that ag is a curve. 
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Figure 1: Mapping a 1-D fiber 

5 Implementing the algorithm 

We now describe the finite element discretization of the algorithm above. Given a 
domain fl C M", we consider a triangulation with interior vertices i^j, j = 1, . . . ,N. 
The nodal functions ip^ are continuous functions which are linear on each element, 
with values on vertices given by ip'^{i'k) — 5jk- The nodal functions span the finite 
element space Vi. 

5.1 Moving Horizontally 

For u'' eViCX = H^{n) and g e L'^{n) cY = H-^{n), we now discretize 

L,{u)z = -Az - P^d2f{-,u)P^z -cQ^z^g- F{u). (6) 



As described in Section |4~T1 this is the main step to identify the fiber ag. 

Functions in X and Y are approximated by elements in Vi, but their identifi- 
cation is different. We take the nodal functions {V'?} as a basis for X^ = Vi^X. 
For u'* £ X'*, we have u^{x) = Y^ j ]i.ji^'j {x) , where u^- = u^{vj). 

For functions g G Pi C 1", we are interested in the values of the (independent) 
functionals ^i(g) = {^p'l,g)o = gi- We take in Vi the dual basis £*,j — l,...,N, 



defined by ii 
coordinates: 



■) = 6ij, so that g{x) — X^iSi ^|(^)' '^^^ mass matrix M changes 



Mu 



M,; 



(V'f,V4)o. 



In coordinates u and g, the expression — Au = g becomes 



Ku: 



K,, = (^f,V'')i, 



where K is the standard stiffness matrix. 

Define inner products (ui,U2)x'' — (Kui,U2) ^^^ (giig2)y'> = (K^^gi,g2) inX^ 
and y '' (here ( , ) denotes the standard inner product in Euclidean space) . Notice 
the isometry (ui,U2)x'' = (gijg2)y'>, where Kuj — gi, i = 1,2. The eigenpairs A^, 
(fi, i £l, have approximations Xf, f^ E X^ obtained by solving 

K^f = AfM<^f, (^f,<£f)x'>=l. 

The approximate eigenfunctions Lp'^ span V^ = V^ C X^ , which may be taken 
arbitrarily close to the vertical subspaces associated to the index set X by choosing 
a small value of h. Similarly, the horizontal subspaces Wx and WV are approximated 
by Wx and WV , the orthogonal complements of V^ and V^ in X and Y respectively. 
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Proposition [5] implies a certain kind of stability. The uniform steepness of fibers 
and the uniform flatness of sheets ensure preservation of flatness. More precisely, if 
F : X = Wx ®Vx ^>-Y = Wy ® Vy is fiat, then it is also flat with respect to the 
decompositions 

X ^Wx® V^, Y = Wy® v^, 

provided that the respective subspaces are sufhciently close to each other. 
Define F^ : X^ ^ Y^ by 

F^{]x) = Ku-M/(u) = F, 

where /(u) is the vector whose coordinates are f{i'j,]i.j)- For small h, F'^ is flat 
with respect to the decompositions X'^ — Wx © V^ and Y'^ — Wy © V^, where 
the horizontal spaces Wx and Wy are orthogonal to V^ — V^ in the discrete inner 
products. 

Assuming z e X'^ in equation ([6]) and taking the L^ inner product of with the 
nodal function tp'^ £ X, we obtain 

(Kz), - {PYd2f{;u)PxZ,^'^)a - c{QxZ,^j'l)„ = g, -F,. 

Since d2f{;u)PxZ G L2(^)^ (P^a2/(-, u)Px^>?)o = (52/(-, u)Fx^, PxV-f )o and we 
are left with discretizing Px ~ I — Qx- In coordinates, Q^z = /_^(z, ^'l) x>^ v'l- 

k 

Once we write z — '^jZjtp'j, the discretization L'' of Lc{u) is expressed in terms 
of the inner products 

{^^,ip1)u {d2fi;u)i;'l,i;i)„, {d^fi-^u),^!;,^^),, {d2fi;u)^1,^^,l, 

where j^k — 1, ■ ■ ■ , N and i, i' G I. In our computations, we replaced 92/(-, u) by 
the vector with coordinates f{vj,'a.j). 

The discretization of the the uptdating m„ i— > m„+i defined in ([5]) becomes then 

u:=u + Px?7, where L*" ?] == ,g - i^''(u). 

5.2 Moving along a fiber 

The finite dimensional inversion of a computable function, as J-g — Qy o F o Tig, is 
not a trivial issue. What is needed is a solver which takes into account the special 
features of maps between vertical subspaces (or, more geometrically, from fibers to 
vertical subspaces). In the examples below, except for the last one, |I[ = f . 

The fact that there was a finite dimensional reduction for the equation F{u) = g 
was implicit in [3], restated in [TB] and stated in a very explicit form (Theorem 2.1) 
in [TB] . Smiley and Chun ( [5D] ) considered the numerical inversion of restrictions of 
F to given fibers, using an inversion algorithm they developed for locally Lipschitz 
maps between Euclidean spaces ([H]). 

As usual, the more we know about F, the sturdier the numerics. The Ambrosetti- 
Prodi case is rather simple: the nonlinearity / interacts only with Ai and 82 f > 0. 
The map F sends fibers to folded vertical lines: as the height of a point in the 
fiber goes from —00 to 00, the height of its image goes monotonically from —00 
to a maximal point and then decreases monotonically to —00. Dropping convexity 
allows for loss of monotonicity, but not of asymptotic behavior, as we shall see in 
the examples of the next section. 

There are theoretical results ([H]) that guarantee that under different, but strin- 
gent, hypotheses the Ambrosetti-Prodi pattern along fibers carries through. The 
numerics may be performed in more general conditions, providing strong evidence 
to the eventual outcome. 
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6 Numerical Examples 

All the examples in this chapter relate to the autonomous equation 



with Dirichlet boundary conditions on the rectangle fl 
smallest three (simple) eigenvalues are 



[0,1] X [0,2], for which the 



^1 = 4 



12.34, A2 



27r^ 



13 

19.74, A3 = —TT 



32.07. 



The nonlinearities / are always appropriate functions. When / is convex, we take 
f'{x) = aarctan(x) +/? for different choices of the asymptotic parameters a and /3. 

Recall that first, given a horizontal affine subspace v + Wx and a right hand 
side g E Y, the algorithm searches for a point Ug E v + Wx in the fiber ag, using 
the iteration described in Section HTTl In each step we solve Equation (O: in the 
examples below, c = 0. Then inversion oi F : ag ^i- g + Vy with basepoint Ug 
obtains, in principle, all solutions of the equation. 

The triangulation was generated with Matlab's PDE Toolbox and the matrices 
were programmed from scratch and compared to those computed by the toolbox, 
whenever possible. 

6.1 Finding Ug in ag 

Consider the Ambrosetti-Prodi situation with f'{x) = aarctan(a;) + /3 satisfying 
lim^j:^±oo f (x) = Ai±(A2-Ai)/2. The right-hand side g(a;) = -100x{x-l)y{y-2) 
is chosen to resemble a very negative multiple of ipi. 

Usually one or two iterations of the horizontal step lead to an error which can 
only decrease by choosing a finer triangulation. Newton's iteration was very suc- 
cessful: continuation arguments were not necessary. An 771-triangulation Tm splits 
each interval [0, 1] and [0, 2] in 2™ equal subintervals. For uq = 100(y92, we present 
the normalized horizontal errors e„ — WPvig — -P"('"n))||/||-fV(<7 ~ ^("o))l|i n = 1,2 



and 3, for triangulations with m = 3, 4 and 5 for the H and H norms. 



m 


ei{H-') 


62 


63 


3 


1.42E-2 


5.27E-5 


4.48E-8 


4 


1.70E-2 


1.12E-4 


3.93E-8 


5 


1.75E-2 


1.31E-4 


4.25E-8 



m 


ei (ff") 


62 


63 


3 


1.97E-2 


9.37E-5 


7.45E-8 


4 


2.36E-2 


1.74E-4 


1.21E-7 


5 


2.44E-2 


1.93E-4 


l.llE-7 



In Figure [2] we show g and the function U3. 





Figure 2: A right-hand side g and Ug e ag obtained from uq = 0. 
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6.2 Finding Solutions: Moving Along a Fiber 

In this section, we prescribe a nonlinearity /, a point uq and study the restriction 
of F to the fiber ag for g — F{uq). The eigenfunctions ipk are normalized in the 
iJ-'^-norm (resp. H^^) in the domain (resp. counter-domain). 

Each example starts with two graphs. In the first, we plot /' and mark with 
dotted lines the relevant eigenvalues. The second graph plots the height of F{u) 
against the height of a point u ^ ag. Informally, it shows how the image of a fiber 
goes up and down: in particular, it indicates the number of solutions of F{u) — 
F{uq). Additional solutions to the equation are then presented. 

6.2.1 Dolph-Hammerstein 
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Figure 3: /' strictly below Ai 
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Figure 4: /' between Ai and A2 

The left of Figure[3]is the graph of /': it lies below the first eigenvalue. The first 
three eigenvalues are marked as dotted lines. The graph on the right illustrates the 
fact that as we move up along the fiber ui = uq + ^^\i the corresponding point in 
the range Fi = F(u\) also moves up. 

Similarly, in Figure IH the derivative of / lies strictly between Ai and A2. Here, 
moving up in the fiber, corresponds to moving down in the range. The graphs are 
consistent with the fact that i^ : X — >■ F is a global diffeomorphism in both cases. 

6.2.2 Ambrosetti-Prodi 

We now return to the example of Section [6. 1[ in which Ai is the only eigenvalue in 
/'(M). Naively, Figures [3] and |4] indicate that, as we move up in the fiber, the image 
under F initially goes up, then down, as shown in Figure O 
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Figure 5: /'(M) n cr(-A) = {AJ = I 

Again, the picture is in agreement with the Ambrosetti-Prodi theorem: below a 
certain height, a point in the vertical line through F{uq) has two preimages. The 
two preimages of F{uq) are shown in Figure [6l 













Figure 6: Ambrosetti-Prodi solutions 



6.2.3 Non-convex /, X = {1} 
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Figure 7: Non-convex / 

Things get more interesting if we relax the condition that / be convex. In 
Figure [7] we analyze the situation in which a non-convex /' interacts only with Ai. 
For uo = —50(^1 + 10(^2 and g = F(uo), the equation F(u) = g has three distinct 
solutions, displayed in Figure [H 
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Figure 8: The three solutions 

The frames in Figure [9] show that the action of F on fibers is not homogeneous. 
The plots show the images under F of fibers ag- with gi = F{—50lpi + Ci(p2), for 
d = 10 (same as Fig. [71), C2 = 45 and C3 — 100. 
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Figure 9: Fibers getting mapped non-uniformly 



6.2.4 Convex /, X = {2} 
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Figure 10: /'(M) n cr(-A) = {A2} = I 



We take / convex, Ran /' = (A2 



A2-A 



^ A2 + 



A2-A 



i) and Mo = —50ip2 + 10ipi. 



2 ' "^ ' 2 
In Figure [TOl heights along the fiber and its image are measured with respect to the 
second eigenfunction ip2- Now, for g = F(uo), there are three preimages, shown in 
Figure [TT] Numerical evidence suggests uniform action of F across fibers. 



6.3 A two dimensional fiber: X = {1, 2} 

We now try to visualize the action of _F on a two-dimensional fiber. 
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Figure 11: Convex f, I — {2}: the three solutions 
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Figure 12: Convex /, X = {1, 2} 



More specificaUy, we examine the fiber ao through the zero function, for which 
F{0) = 0. Consider the circle C in Vx, the vertical plane spanned by ipi and (/?2, 
shown in Figure [T51 Let Ca C ao be the curve HqCC), which projects bijectively 
under Qx to C. 
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Figure 13: Solutions U and D on the circle C and images. 

The fish-shaped curve in Figure [T^ is the projection of F{Ca) under Qy in V^. 
Seven points and their images were given common labels. Let g, marked with a 
bullet, be the point of self-intersection of this curve. Clearly g has two preimages U 
and D between points 2 and 3 and 6 and 7, respectively. Radial lines in the domain 
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from the origin to points in C give rise to lines from -F(O) == to points in F(Ca), 
as seen in Figure 1141 We then obtain two approximate preimages L and R along 
the horizontal axis. 

The four approximate preimages were then taken as initial guesses for Newton's 
Method and the four computed solutions are illustrated in Figure [TSl 



150 






; n_, 








60 
















' 


' 


' 


' 


100 








■^3 
"2 








40 










50 



-50 


w w 
■■- • - - 


w^ w 





R 








20 

Li_ 

-20 






™4 ^3 




L 




=, 




S 


e^ 


''4 


S 


^2 z 


s. 








Sj 




















100 






s 










-40 
















s . 








-60 












-150 


-100 


-50 
F1 





-1 


)0 -100 


-50 


50 


100 


1f 















u 


1 



















Figure 14: Solutions L and R on the ul axis 









(a) U and D 












(b) L and R 



Figure 15: Computed Solutions, 2-D case 
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